Kakava et al scLepr hypothalamus paper figures
**What this notebook does**
This notebook reproduces the figures of the analysis the tdtomato sorted Lepr+ cells from the hypothalamus by Kakava at al
**Where does the data come from?**
The shaping, preprocessing and the analysis of the single cell data is provided in seperate files. Here, h5ad files containing this information is loaded and used to reproduce the figures in the manuscript. The raw data including the fastq files and count tables can be downloaded the data from GEO (GSE65970, will soon be publicly avaialable).
We provide the raw count table that can be used to regenerate all the figures of hte manuscipt here on Guthub. Following the notebooks, you can preprocess the data to generate the anndata objects, saved as h5ad files, that are used for this study under the './rawdata' folder. These files are too large to load
Pyscenic analysis are performed on the high performance cluster of hte Utrecht Bioinformatics center
**Notes**
The downsampled dataset is used for the majority of the analysis, except for calling the active gene regulatory networks. This dataset is normalised, but not regressed for nCounts
'tomcells' refer to the whole dataset, where only tdToamto+ sorted cells are used. Neurons are subsetted using 'neuronal cluster' information based on Tubb3 and Mapt expression. Nonneuronal cells are analysed seperately.
Note that the cell type names are slightly different. 'Chrh' neurons are naes to 'Trh' neurons. Similarly, the order of endothelial cluster numbers are changed in the manuscript. These are minor changes, however, we will update to the latest version upon publication of the manuscript
In [115]:
import numpy as np
import pandas as pd
import scanpy as sc
import os
sc.settings.verbosity = 0 # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_versions()
data = './'
In [116]:
sc.settings.set_figure_params(dpi=80)
upload the data¶
Upload the already analysed dataset
In [117]:
adata = sc.read('./rawdata/Lepr_hypothalamus_tomcells_analysis.h5ad')
adata_down = sc.read('./rawdata/Lepr_hypothalamus_tomcells_down_analysis_celltype.h5ad')
In [118]:
adata_down
Out[118]:
some generic stats¶
Here, we provide some of the statistics used in the manuscript
Some statistics for cell numbers:¶
In [119]:
for cell_type in adata_down.obs['cell_types'].unique():
print(cell_type,": ",len(adata_down[adata_down.obs['cell_types'].isin([cell_type])].obs_names))
The difference between the ERCC content of fed and fasted conditions of Agrp neurons¶
In [120]:
adata_fasted = adata[adata.obs['condition'] == 'fasted']
adata_fed =adata[adata.obs['condition'] == 'fed']
In [121]:
print('ERCC_genes fasted',adata_fasted[adata_fasted.obs['cell_categories'] == 'Agrp_neurons'].obs['ERCC_genes'].median())
print('ERCC_genes fed',adata_fed.obs[adata_fed.obs['cell_categories'] == 'Agrp_neurons']['ERCC_genes'].median())
print('n_counts fasted',adata_fasted[adata_fasted.obs['cell_categories'] == 'Agrp_neurons'].obs['n_counts'].median())
print('n_counts fed',adata_fed[adata_fed.obs['cell_categories'] == 'Agrp_neurons'].obs['n_counts'].median())
In [122]:
sc.pl.violin(adata[adata.obs['cell_categories'] == 'Agrp_neurons'],
keys = ['ERCC_genes'], groupby='condition')
sc.pl.violin(adata[adata.obs['cell_categories'] == 'Agrp_neurons'],
keys = ['n_counts'], groupby='condition')
Distribution of n_counts, conditions and clusters in original and downsampled datasets¶
In [123]:
sc.pl.umap(adata, color=['condition','n_counts'], use_raw=True,save="_tomcells_condition_counts.pdf")
sc.pl.umap(adata_down, color=['condition','n_counts'], use_raw=True,save="_tomcells_condition_counts.pdf")
The distribution of clusters
In [124]:
sc.pl.umap(adata, color=['leiden'], use_raw=True,save="_tomcells_leiden.pdf")
sc.pl.umap(adata_down, color=['leiden'], use_raw=True,save="_tomcells_down_leiden.pdf")
Figure 1¶
Fig. 1b : Cell types¶
In [125]:
sc.pl.umap(adata_down, color=['cell_types'], use_raw=True,
save="_tomcells_celltypes_Fig1b.pdf")
Fig. 1c : Cell types violin plots¶
Violin plot showing selected lineage marker expression by Leiden clusters¶
In [126]:
ax = sc.pl.stacked_violin(
adata_down, ['Tubb3','Mapt','Rbfox3','Slc17a6','Slc32a1',
'Pecam1','Agt','Gfap','Rax','Col1a2','Lepr','Olig2','Plp1',
'Cd53','Tshb','Nrn1'], groupby='cell_types', swap_axes=True,
var_group_labels=['Tubb3','Mapt'], dendrogram=True,
save='_tomcells_celltypes_Fig1c_corrected.pdf')
Supp fig S1b¶
Batches and conditions used (Related to Figure 1)
In [127]:
sc.pl.umap(adata_down, color=['plate', 'condition'],
use_raw=True,save="_batches_condition_FigS1b.pdf")
Supp fig S1c¶
Quality control measures (Related to Figure 1)
In [128]:
adata
Out[128]:
In [129]:
sc.pl.violin(adata_down, keys=['n_genes','n_counts','percent_mito','ERCC_genes'],
jitter=0.4, multi_panel=True,save="_tomcells_qcmetrics_figS1c.pdf")
sc.pl.scatter(adata_down, x='n_counts', y='n_genes',color = 'percent_mito',
save= "_tomcells_complexity_figS1c.pdf")
Supp fig 2a¶
Hierarchical clustering heatmap of cell types (Related to Figure 1)
In [130]:
#Clustermap computed by sc.tl.dendogram
sc.pl.correlation_matrix(adata_down,groupby='leiden',figsize=[10,10],dendrogram ='dendogram_leiden',
save='_tomcells_down_dendogram_leiden_figS2a.pdf')
Supp fig 2b¶
Lepr expression (Related to Figure 1)
In [131]:
sc.pl.umap(adata_down, color=['Lepr'], use_raw=True,save="_tomcells_sup_Lepr.pdf")
Figure 2¶
Upload the data of neuronal and nonneuronal populations
In [132]:
adata_down_neuron = sc.read('./rawdata/Lepr_hypothalamus_tomcells_analysis_neuron.h5ad')
adata_down_nonneuron = sc.read('./rawdata/Lepr_hypothalamus_tomcells_down_analysis_nonneuron.h5ad')
In [133]:
adata_down_neuron
Out[133]:
In [134]:
adata_down_nonneuron # Excludes pitutiary
Out[134]:
Fig. 2a¶
Neuronal cell types displayed on UMAP
In [135]:
sc.pl.umap(adata_down_neuron, color=['condition'])
In [136]:
sc.pl.umap(adata_down_neuron, color=['cell_types']
, use_raw=True,save="_tomcells_neuron_celltypes_fig2a.pdf")
Fig. 2b¶
Neuronal marker gene expression by neuronal cell types
In [137]:
list_neuronal = ['Slc32a1', 'Slc17a6','Agrp', 'Npy', 'Pomc', 'Cartpt', 'Kiss1', 'Pdyn',
'Tac2','Ghrh', 'Gal', 'Trh', 'Th', 'Htr2c', 'Esr1', 'Ghr', 'Crhr1', 'Irs4',
'Ghsr', 'Egr1', 'Sst', 'Prlr', 'Pgr', 'Gfra1', 'Reln']
list_neuronal= np.intersect1d(list_neuronal,adata_down_neuron.raw.var_names)
list_neuronal
Out[137]:
Stacked violin plots. Markers are displayed alphabetically; 'Slc32a1' and 'Slc17a6' manually moved to the top using the Adobe Illustrator
In [138]:
ax = sc.pl.stacked_violin(
adata_down_neuron, list_neuronal, groupby='cell_types', swap_axes=True,
var_group_labels=['0'], dendrogram=True,
save='_neuronal_down_Fig2b_list_neuronal_genes_corrected.pdf')
Fig. 2c¶
In [139]:
sc.pl.umap(adata_down_neuron, color=['Slc32a1', 'Slc17a6'],
use_raw=True,save="_neurons_Fig2c.pdf")
Fig. 2d¶
In [140]:
sc.pl.umap(adata_down_neuron, color=['Th', 'Ddc'],
use_raw=True,save="_neurons_Fig2d.pdf")
Fig. 2e¶
In [141]:
sc.pl.umap(adata_down_neuron, color=['Crhr1', 'Crhr2'],
use_raw=True,save="_neurons_Fig2e.pdf")
Fig. 2f¶
In [142]:
sc.pl.umap(adata_down_neuron, color=['Npy1r', 'Npy2r'],
use_raw=True,save="_neurons_Fig2f.pdf")
Supp Fig. S3a¶
Neuropeptide expression by Lepr+ neurons (related to fig 2)
In [143]:
genes = ['Adcyap1','Adm','Agrp','Apln','Avp','Calca','Cartpt','Cbln1',
'Cbln2','Cbln3','Cbln4','Cck','Chga','Chgb','Cort','Crh',
'Dbi','Edn1','Edn2','Edn3','Gal','Galp','Gast','Ghrh',
'Gnrh1','Grp','Hcrt','Igf1','Igf2','Insl6','Kiss1','Nms',
'Nmu','Npvf','Npy','Nts','Oxt','Oxt','Pdyn','Penk',
'Pmch','Pnoc','Pomc','Prlh','Prok2','Pthlh','Retn','Rln1',
'Tac1','Tac2','Trh','Uts2','Vgf','Vip']
genes= np.intersect1d(genes,adata_down_neuron.raw.var_names)
genes
Out[143]:
In [144]:
ax = sc.pl.stacked_violin(
adata_down_neuron, genes, groupby='cell_types', swap_axes=True,
dendrogram=False,
save='_neurons_FigS3a_neuropeptides.pdf')
Supp Fig. S3b¶
Neuropeptide receptor expression by Lepr+ neurons (related to fig 2)
In [145]:
receptors = ['Adcyap1r1','Avpr1a','Avpr1b','Avpr2','Cckar','Cckbr','Crhr1','Crhr2','Esr1','Esr2',
'Galr1','Galr2','Galr3','Ghr','Ghsr','Gnrhr','Grpr','Irs4','Kiss1r','Lepr',
'Mc3r','Mc4r','Mchr1','Npffr1','Npy1r','Npy2r','Npy5r','Ntsr1','Ntsr2',
'Oprd1','Oprk1','Oprl1','Oprm1','Oxtr','Pgr','Prokr1','Prokr2','Rxfp1',
'Rxfp2','Rxfp3','Sstr1','Sstr2','Sstr3','Sstr4','Sstr5','Tacr1','Tacr2',
'Tacr3','Trhr','Vipr1','Vipr2']
receptors= np.intersect1d(receptors,adata_down_neuron.raw.var_names)
receptors
Out[145]:
In [146]:
ax = sc.pl.stacked_violin(
adata_down_neuron, receptors, groupby='cell_types', swap_axes=True,
dendrogram=False,
save='_neurons_FigS3b_neuropeptides_receptors.pdf')
Figure 3¶
Analysis of Trh neuronal clusters
In [147]:
adata_down_neuron.obs['cell_types']
Out[147]:
Fetch Trh neurons
In [148]:
adata_down_trhneurons = adata_down_neuron[adata_down_neuron.
obs['cell_types'].isin(["Crhr2_Npw_neurons",
"Crhr2_Prlh_neurons","Chrh2_Crh_neurons"])].copy()
#.str.contains('Crhr2')] -> Correct misspelled 'Chrh2'
In [149]:
adata_down_trhneurons
Out[149]:
Recalculate the neighborhood, umap embedding
In [150]:
sc.pp.neighbors(adata_down_trhneurons, n_neighbors=15, n_pcs=40) #
sc.tl.umap(adata_down_trhneurons,spread = 1.5,maxiter =1000,min_dist=1.5)
Fig. 3a¶
Marker gene expression on dotplots
Note that the Crhr2 clusters are renamed to Trh
In [151]:
Trh_genes = ['Trh', 'Lhx1', 'Lhx1os', 'Ptger3', 'Adra1b',
'Cxcl12', 'Cnr1', 'Tac1', 'Pirt', 'Crhr2','Prlr']
In [152]:
ax = sc.pl.dotplot(
adata_down_trhneurons, Trh_genes, groupby='cell_types',
dendrogram=False,figsize = [5,2],
save='_trhneurons_Fig3a_left.pdf')
In [153]:
sc.pl.umap(adata_down_trhneurons, color=['cell_types'],
use_raw=True,save="_trhneurons_Fig3a_right.pdf")
# The colors dont match the figure, but this is acceptable
Fig. 3b¶
In [154]:
sc.pl.umap(adata_down_trhneurons, color=['Lhx1', 'Lhx1os', 'Ptger3', 'Adra1b' ],
use_raw=True,save="_trhneurons_Fig3b.pdf")
Fig. 3c¶
In [155]:
sc.pl.umap(adata_down_trhneurons, color=['Cnr1', 'Gal', 'Nts', 'Crh'],
use_raw=True,save="_trhneurons_Fig3c.pdf")
Figure 4¶
Nonneuronal cells
Fig 4a¶
In [156]:
tanycytes = adata_down_nonneuron[adata_down_nonneuron.obs['cell_types'].str.contains('Tanycyte')]
tanycytes # 156 cells
Out[156]:
stacked violin plots of tanycyte markers
In [157]:
tanycyte_markers = ['Rax', 'Vim', 'Ppp1r1b', 'Nes', 'Gfap','Crym', 'Frzb', 'Col25a1', 'Adm',
'Cacna2d2', 'Vcan', 'Cd59a', 'Slc17a8','Tll1', 'Lrp2','Cntfr']
In [158]:
ax = sc.pl.stacked_violin(
tanycytes, tanycyte_markers, groupby='cell_types', swap_axes=True,
dendrogram=False,
save='_tanycytes_Fig4a.pdf')
Fig 4b¶
Umaps of Tanycyte clusters with: Rax, Vim, Ppp1r1b, Nes, Gfap,Crym, Frzb, Col25a1, Adm, Cacna2d2, Vcan, Cd59a, Slc17a8,Tll1, Lrp2,Cntfr.
In [159]:
endothelial = adata_down_nonneuron[adata_down_nonneuron.obs['cell_types'].str.contains('endothelial')]
endothelial # 150 cells
Out[159]:
dotplots of endothelial marker genes
In [160]:
genes = ['Pecam1', 'Ly6c1', 'Ly6a', 'Slco1c1','Slc38a5', 'Vwf', 'Slc16a1', 'Vegfc',
'Sema3g', 'Gkn3', 'Acta2', 'Myh11', 'Pdgfrb', 'Vtn']
In [161]:
ax = sc.pl.dotplot(
endothelial, genes, groupby='cell_types', dot_max=1, color_map = 'coolwarm',
dendrogram=False,
save='_endothelial_Fig4b.pdf')
Fig 4c:¶
Violin plots of VLMC clusters with: Igf2, Gnrh1, Penk, Spp1
In [162]:
VLMC = adata_down_nonneuron[adata_down_nonneuron.obs['cell_types'].str.contains('Vlmc')]
VLMC # 63 cells
Out[162]:
stacked violin plots of VLMC marker genes
In [163]:
genes = ['Igf2', 'Gnrh1', 'Penk', 'Spp1']
In [164]:
ax = sc.pl.stacked_violin(
VLMC, genes, groupby='cell_types', swap_axes=True,
dendrogram=False,figsize= [3,4],
save='_VLMC_Fig4c.pdf')
Fig 4d¶
Violin plots of glial (Oligo, OPC, Astro, MM) clusters marked with: Mobp, Ermn, Mbp, Fyn, Bmp4, Agt, Gfap, Itih3, Cx3cr1, Mrc1.
In [165]:
astroglial = adata_down_nonneuron[adata_down_nonneuron.obs['cell_types'].isin(['Microglia','OPCs',
'oligodendrocytes','Astrocytes'])] #eclude 'Nrn1_non_neuronal' unknown
astroglial # 39 cells
Out[165]:
stacked violin plots of glial markers
In [166]:
genes = ['Mobp', 'Ermn', 'Mbp', 'Fyn', 'Bmp4', 'Agt', 'Gfap', 'Itih3', 'Cx3cr1', 'Mrc1']
In [167]:
ax = sc.pl.stacked_violin(
astroglial, genes, groupby='cell_types', swap_axes=True,
dendrogram=False,figsize=[2,5],
save='_astroglial_Fig4d.pdf')
Fig 4e¶
Umaps showing astrocytes only. Markers: Agt, Gfap, Itih3
In [168]:
astrocytes = adata_down_nonneuron[adata_down_nonneuron.obs['cell_types'].isin(['Astrocytes'])]
astrocytes # 21 cells
Out[168]:
In [169]:
genes = ['Agt', 'Gfap', 'Itih3']
In [170]:
sc.pl.umap(astrocytes, color=genes, size= 5000,
use_raw=True,save='_astroglial_Fig4e.pdf')
Figure S4¶
comparison of cell types to Chen et al
Export the rawdata (norm, log-transformed), hvgenes subsetted data and cell type annotation as csv. Combine them manually
In [171]:
rawdata = pd.DataFrame(adata_down.raw.X,index=adata_down.raw.obs_names, columns=adata_down.raw.var_names)
rawdata
Out[171]:
save the data
In [172]:
rawdata.T.to_csv('./cell_type/tomcells_rawdata_normlog.csv')
rawdata[adata_down.var_names].T.to_csv('./cell_type/tomcells_rawdata_normlog_hvgenes.csv')
save the cell type annotation
In [173]:
cell_types = pd.DataFrame(adata_down.obs['cell_types']).T
cell_types.to_csv('./cell_type/tomcells_cell_types.csv')
cell_types
Out[173]:
SciBet online tool is used to compare the Lepr+ cell identity to the published cell types from Chen et al. The reference dataset is provided on the website as one of the default datasets
In [174]:
predicted_cell_type = pd.read_csv('./rawdata/cell_type/Chen/predicted_cell_type.csv',index_col=0)
probability_matrix = pd.read_csv('./rawdata/cell_type/Chen/probability_matrix.csv',index_col=0)
In [175]:
adata_down.obs['predicted_cell_type'] = predicted_cell_type
In [176]:
probability_matrix = sc.AnnData(probability_matrix)
probability_matrix.obs['cell_types'] = adata_down.obs['cell_types']
probability_matrix
Out[176]:
In [177]:
probability_matrix.var_names
Out[177]:
In [178]:
sc.pl.matrixplot(probability_matrix, groupby='cell_types',var_names=probability_matrix.var_names,
dendrogram=True, save='_cell_types_compared_to_Chen_figS4.pdf')
Out[178]:
Figure 5¶
pyscenic analysis
Load the anndata files with scenic analysis
In [179]:
adata_scenic_down = sc.read('./rawdata/pyscenic/tomcells/tomcells_down_pyscenic.h5ad')
In [180]:
adata_scenic_down
Out[180]:
Fig 5a¶
cell categories displayed on t-SNE space calculated based on regulon activity
In [181]:
sc.pl.tsne(adata_scenic_down, color=['cell_categories'],save='_cell_categories_Fig5a.pdf')
Fig 5b¶
Hieararchical clustering of cell types baed on regulon activity
Load hte data file that contain the information on regulons for tomcells, calculated using log normalised data to maintain as much information as possible
In [182]:
adata_auc_mtx = sc.read('./rawdata/pyscenic/tomcells/tomcells_auc_mtx.h5ad')
In [183]:
adata_auc_mtx
Out[183]:
Plot the dendrogram. Note that the differentially active regulons were added using Adobe Illustrator
In [184]:
sc.pl.dendrogram(adata_auc_mtx,groupby = 'cell_types',save='_tomcells_auc_mtx_dendrogram_Fig5b.pdf')
Out[184]:
For instance, the differences between neural and nonneural clusters are calcualted as follows:
In [185]:
sc.tl.rank_genes_groups(adata_auc_mtx,groupby='level_1',
groups=['nonneural'],reference ='neural',method='wilcoxon') #wilcoxon
sc.pl.rank_genes_groups(adata_auc_mtx,groupby='level_1',groups=['nonneural'],n_genes=40)
Fig S5a¶
conditions displayed on t-SNE space calculated based on regulon activity
In [186]:
sc.pl.tsne(adata_scenic_down, color=['condition'],save='_condition_FigS5a.pdf')
Fig S5b¶
dotplot displaying regulon activity in cell types
In [187]:
sc.pl.dotplot(adata_auc_mtx,var_names=adata_auc_mtx.var_names,groupby='cell_types',
color_map='Spectral',dendrogram=True, save='_tomdata_regulon_activity_figS5b.pdf')
Out[187]:
Figure 6¶
Comparison of fed and fasted Agrp neurons
Load the Agrp neurons data, that was reanalysed from scratch. The raw data was downsampled to the median number of reads
In [188]:
agrpneurons = sc.read('./rawdata/Lepr_hypothalamus_agrpneurons_down_preprocess.h5ad',
compression='gzip')
In [189]:
agrpneurons.obs['n_counts'].median()
Out[189]:
Following downsampling to the median, cell cluster largely together. However, cells from the fasted condition are not equally distributed
In [190]:
sc.pl.umap(agrpneurons,color = 'condition',add_outline=True,
save = '_agrpneurons_down_condition_lined.pdf')
First, we will look at the data a bit closer. Some of these stats are used in the paper and/or provided as supplementary tables
Calculate the differentially expressed genes between conditions using wilcoxon rank-sum test
In [191]:
sc.tl.rank_genes_groups(agrpneurons,'condition', method='wilcoxon',n_genes=500)
Generate a table with logfoldchanges, scores and corrected p-calues for top500 differentially expressed genes
In [192]:
result = agrpneurons.uns['rank_genes_groups']
groups = result['names'].dtype.names
result_DGE = pd.DataFrame(
{group + '_' + key[:1]: result[key][group]
for group in groups for key in ['names', 'scores','logfoldchanges','pvals_adj']})
result_DGE
Out[192]:
Genes upregulated upon fasting with a FDR <0.01
In [193]:
result_DGE[result_DGE['fasted_p'] < 0.01] #263 genes FDR<0.01
Out[193]:
Genes downregulated upon fasting with a FDR <0.01
In [194]:
result_DGE[result_DGE['fed_p'] < 0.01] #35 genes FDR<0.01
Out[194]:
Fig 6a¶
The following are selected from the fasted-enriched genes
In [195]:
#order manually. Keep 'Lepr', 'Agrp', 'Npy' on the top, the rest follows the p-value (FDR)
agrp_fasting_genes=['Lepr', 'Agrp', 'Npy', 'Vgf','Ghr','Acvr1c','Mt1','Trpm3','Asb4',
'Sorcs1', 'Slc8a1','Fkbp5','Sema3a','Slc4a4','Ghsr','Trpv2','Arrb1']
In [196]:
ax = sc.pl.dotplot(
agrpneurons, agrp_fasting_genes, groupby='condition',
dendrogram=False, color_map= 'coolwarm', #'Spectral_r',#figsize = [5,2],
save='_agrp_fasting_genes_Fig6a_coolwarm.pdf')
Fig 6b¶
The folowing are the top20 genes enriched in fed cells
In [197]:
fed_genes = result_DGE[result_DGE['fed_p'] < 0.01]['fed_n'][:20]
In [198]:
sc.pl.dotplot(
agrpneurons, fed_genes, groupby='condition',
dendrogram=False, color_map= 'coolwarm', #'Spectral_r',#figsize = [5,2],
save='_agrp_fed_genes_Fig6b_coolwarm.pdf')
Out[198]:
Fig 6c¶
Regulon enrichment in cells from fasted versus fed conditions. For agrp neurons, we have recalculated the regulon activity enrichment for single cells using the data from agrp neurons which was downsampled to the median
The regulons detected here are not the same as those detected in the whole dataset. Part of this maybe because only Agrp neurons are analysis. Regulons enriched in every cell is calculated in comparison to other cells, and we may see more active regulons if this dataset is analysed in comparison to all Lepr cells or the whole hypothalamus. Here, we focus on the differences between conditions or subclusters
Load the data
In [199]:
agrpneurons_pyscenic = sc.read('./rawdata/agrpneurons_pyscenic.h5ad')
In [200]:
agrpneurons_pyscenic
Out[200]:
Load the active regulons as a data frame
In [201]:
auc_mtx_down = pd.read_csv('./rawdata/pyscenic/agrpneurons/agrpneurons_down_to_median.auc.csv', index_col=0)
In [202]:
auc_mtx_down
Out[202]:
convert into anndata
In [203]:
adata_auc_mtx_down = sc.AnnData(auc_mtx_down)
adata_auc_mtx_down.obs['condition'] = adata_down.obs['condition']
In [204]:
adata_auc_mtx_down
Out[204]:
In [205]:
sc.tl.rank_genes_groups(adata_auc_mtx_down,groupby = 'condition',method='wilcoxon')
sc.pl.rank_genes_groups(adata_auc_mtx_down,groups = ['fasted'],fontsize = 15,
save='_fasted_fed_regulons_wilcoxon_Fig6c.pdf')
Figure S6a¶
Gene regulatory network aanalysis of Agrp neurons
In [206]:
import seaborn as sns
import matplotlib as mpl
import matplotlib.pyplot as plt
from pyscenic.rss import regulon_specificity_scores
Generate scaled dataframe for regulon activity¶
In [207]:
auc_mtx_Z_down = pd.DataFrame( index=auc_mtx_down.index )
for col in list(auc_mtx_down.columns):
auc_mtx_Z_down[ col ] = ( auc_mtx_down[col] - auc_mtx_down[col].mean()) / auc_mtx_down[col].std(ddof=0)
Generate a heatmap¶
In [208]:
# From the pyscenic analysis pipeline
def palplot(pal, names, colors=None, size=1):
n = len(pal)
f, ax = plt.subplots(1, 1, figsize=(n * size, size))
ax.imshow(np.arange(n).reshape(1, n),
cmap=mpl.colors.ListedColormap(list(pal)),
interpolation="nearest", aspect="auto")
ax.set_xticks(np.arange(n) - .5)
ax.set_yticks([-.5, .5])
ax.set_xticklabels([])
ax.set_yticklabels([])
colors = n * ['k'] if colors is None else colors
for idx, (name, color) in enumerate(zip(names, colors)):
ax.text(0.0+idx, 0.0, name, color=color, horizontalalignment='center', verticalalignment='center')
return f
In [209]:
cats = sorted(list(set(adata_down.obs['condition'])))
In [210]:
colors = sns.color_palette('bright',n_colors=len(cats) )
colorsd = dict( zip( cats, colors ))
colormap = [ colorsd[x] for x in adata_down.obs['condition'] ] #cellAnnot['Celltype_Garnett']
Generate a color map
In [211]:
sns.set()
sns.set(font_scale=0.8)
fig = palplot( colors, cats, size=1.0)
Calculate regulon specificity scores based on condition
In [212]:
rss_condition_down = regulon_specificity_scores( auc_mtx_down, agrpneurons.obs['condition'] )
rss_condition_down
Out[212]:
Calculate the rop regulated regulons using the regulon specificity scores
In [213]:
topreg_down = []
for i,c in enumerate(cats):
topreg_down.extend(
list(rss_condition_down.T[c].sort_values(ascending=False)[:5].index)
)
topreg_down = list(set(topreg_down))
In [214]:
sns.set()
sns.set(font_scale=1.0)
sns.set_style("ticks", {"xtick.minor.size": 1, "ytick.minor.size": 0.1})
g = sns.clustermap(auc_mtx_Z_down[topreg_down], square=False, linecolor='gray',
yticklabels=True, vmin=-2, vmax=6, row_colors=colormap,
cmap="YlGnBu", figsize=(21,16))
g.ax_heatmap.set_ylabel('Cells')
g.ax_heatmap.set_xlabel('Regulons')
g.ax_row_colors.set_yticks([0.5])
g.ax_row_colors.set_yticklabels(['condition'])
g.cax.set_visible(False)
g.savefig("./figures/heatmap_regulons_agrpneurons_down_topreg_rowclust_figS6a.pdf")
Fig. S6b¶
In [215]:
sc.pl.stacked_violin(adata_auc_mtx_down,var_names=['Rfx3','Thoc2','Hivep3','Foxn3','E2f6',
'Npdc1','E2f6','Jun','Fosb','Hivep3','Hlf'],groupby='condition',use_raw=False,swap_axes=True,
save='_fasted_fed_regulons_fig6b_left.pdf')
sc.pl.stacked_violin(agrpneurons,var_names=['Rfx3','Thoc2','Hivep3','Foxn3','E2f6',
'Npdc1','E2f6','Jun','Fosb','Hivep3','Hlf'],groupby='condition',swap_axes=True,
save='_fasted_fed_TFs_expression_fig6b_right.pdf')
Out[215]:
Figure 7¶
CELLECT analyis on BMI GWAS, to be contributed by the Pers group
Tables¶
In [216]:
result = agrpneurons.uns['rank_genes_groups']
groups = result['names'].dtype.names
cluster = agrpneurons
test = pd.DataFrame(
{group + '_' + key[:1]: result[key][group]
for group in groups for key in ['names', 'scores','logfoldchanges','pvals_adj']})
gene_ids = cluster.raw.var.index.values
#recalculate the averages and pct
obs = cluster.raw[:,gene_ids].X
obs = pd.DataFrame(obs,columns=gene_ids,index=cluster.obs['condition'])
average_obs = obs.groupby(level=0).mean()
obs_bool = obs.astype(bool)
fraction_obs = obs_bool.groupby(level=0).sum()/obs_bool.groupby(level=0).count()
df0 = pd.concat([average_obs.T.add_suffix("_mean"), fraction_obs.T.add_suffix("_pct")],
axis=1).sort_index(axis=1,ascending=True)
df1 = df0.loc[test.iloc[:,0]] #enriched in fasted
df2 = df0.loc[test.iloc[:,4]] #enriched in fed
Get the data and names of the differentially expressed genes in fasted cells
In [217]:
diff_fasted = test[test['fasted_l'] >1.5][test['fasted_p'] < 0.01][test.columns[0:3]]
#diff_fasted = test[test['fasted_l'] >1.5][test['fasted_p'] < 0.01]
diff_fasted.index = range(93)
diff_fasted
Out[217]:
Get the data and names of the differentially expressed genes in fed cells
In [218]:
diff_fed = test[test['fed_l'] >1.5][test['fed_p'] < 0.01][test.columns[4:7]]
#diff_fasted = test[test['fasted_l'] >1.5][test['fasted_p'] < 0.01]
diff_fed.index = range(13)
diff_fed
Out[218]:
In [ ]: